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Abstract: We study how to numerically simulate quantum fermions out of thermal 
equilibrium, in the context of electroweak baryogenesis. We find that by combining the 
lattice implementation of Aarts and Smit [1] with the "low cost" fermions of Borsanyi 
and Hindmarsh [2], we are able to describe the dynamics of a classical bosonic system 
coupled to quantum fermions, that correctly reproduces anomalous baryon number 
violation. To demonstrate the method, we apply it to the 1+1 dimensional axial U{1) 
model, and perform simulations of a fast symmetry breaking transition. Compared to 
solving all the quantum mode equations as in [1] , we find that this statistical approach 
may lead to a significant gain in computational time, when applied to 3+1 dimensional 
physics. 
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1. Introduction 

In electroweak baryogenesis [3] the source of baryon number non-conservation is the 
quantum anomaly of fermions chirally coupled to the Standard Model SU(2) gauge field. 
When the gauge field evolves in such a way that its Chern-Simons number changes, the 
fermion, and hence B(aryon) and L(epton), number changes as 



B{t) - B{0) = L{t) - L(0) = nf[N^{t) - N^], 



(1.1) 



where n/ = 3 is the number of fermion generations in the Standard ModeL The 
question of successful baryogenesis thus reduces to whether a permanent change of 
Chern-Simons number can take place in the early Universe, presumably under the 
influence of CP- violation and the back-reaction of the fermions. 

Various models of baryogenesis have been proposed, of which the most popular 
(and most developed) is "hot" clcctrowcak baryogenesis [4], where walls of bubbles 
nucleated in a first order phase transition interact in a CP-violating manner with the 
fermions in the hot plasma. In this way a net left-right fermion asymmetry is generated 
inside and outside the bubbles, and equilibrium gauge dynamics (sphaleron transitions) 
convert this asymmetry into a baryon asymmetry. 

The rate of sphaleron transitions can reliably be calculated in thermal equilibrium 
using sophisticated Monte-Carlo methods [5, 6]. In such a setup, fermions can be 
included in terms of effective couplings for the bosonic theory, for instance through 
dimensional reduction [7]. 

An alternative scenario is "Cold" electroweak baryogenesis [8, 9, 10, 11, 12, 13], 
where the electroweak phase transition does not involve bubble nucleation, but instead 
a fast quench of the Higgs potential. Here, baryon number violating processes arc not 
equilibrium Sphaleron transitions, but complicated out-of-equilibrium field dynamics. 

Numerical real-time simulations of electroweak baryogenesis have until now ne- 
glected dynamical fermions. Instead, purely bosonic systems are evolved and baryon 
number has simply been assumed to follow the gauge field Chern-Simons number in 
accordance with the anomaly equation, ignoring fermionic backreaction. 

One case where this is certainly not allowed is for minimal electroweak baryogenesis, 
since CP-violation in the Standard Model originates from the fermion mass matrix. A 
possible approach employed in [14, 8, 12] is to integrate out the fermions in the path 
integral or in perturbation theory, thus recovering CP-violation effects in terms of a 
series of higher-dimensional bosonic terms. 

The current understanding that Standard Model CP-violation is strongly sup- 
pressed at high temperatures, and therefore insufficient for successful baryogenesis 
follows from such a computation (see for instance [14, 15]). In contrast, at low tem- 
peratures relevant for "Cold" baryogenesis, recent calculations have shown that the 
suppression is absent [16, 17, 18, 19], and direct numerical simulations have in turn in- 
dicated that Standard Model CP- violation may in fact be large enough to accommodate 
the observed asymmetry [13, 20]. 

A possible caveat to this procedure is that it is based on a gradient expansion 
in the gauge and Higgs fields, which may not be valid during electroweak symmetry 
breaking. And so although the work in [13, 20] is very encouraging indeed, it would 
be even better not having to integrate out the fermions, but include them directly in 



- 2 - 



real-time simulations of the transition. In this way, the CKM matrix and CP- violation 
could be included from first principles. 

In [1] Aarts and Smit showed how to implement quantum fermions in real-time, 
coupled to classical bosonic gauge and scalar fields. The method involves a proper 
lattice discretization in Minkowski space, and the realisation that since fermions are 
bilinear in the action, the field operators can be expanded into mode functions, in terms 
of time-independent creation-annihilation operators. These mode functions can then 
be solved in the classical bosonic background, with the back-reaction on the bosonic 
fields defined as the quantum averages over the creation-annihilation operators for some 
given initial state. 

In practice, the problem is that for every momentum mode k (equal to the number 
of lattice sites , where D is the number of spatial dimensions), one needs to solve a 
separate real-time field equation (the mode function equation) for which the numerical 
effort is also proportional to the number of lattice sites. Hence the total numerical 
problem scales as n^^, and quickly becomes unmanageable for large three-dimensional 
lattices. Large lattices are often required in baryogenesis simulations to accommodate 
extended objects such as sphalerons and for having enough infrared modes for a fast 
quench to be correctly reproduced. 

Some time ago [2], Borsanyi and Hindmarsh showed how to replace the mode 
equations by an ensemble of fermion field realisations, approximating the quantum 
fermion expectation values through a statistical averaging procedure. In the context 
of a scalar-fermion theory, they showed that one can significantly reduce the numerical 
effort, at least in three dimensions. This is because the number of random realisations 
in the ensemble Ng can be much smaller than n^. 

In this work, we will implement the "low cost fermion" or "fermion ensemble" 
method of Borsanyi and Hindmarsh to the 1-1-1 dimensional axial-U(l) model with 
fermions of Aarts and Smit. This will act as a toy model for the electroweak part of 
the Standard Model, and will provide a testing ground for the method. In particular, 
we will investigate whether this method correctly reproduces the anomaly equation, 
charge conservation and the correct dynamics, and determine how large the fermion 
ensemble needs to be to get reliable results. We also want to understand when it is 
correct to neglect fermion backreaction for the boson dynamics. 

The paper is structured as follows: In section 2, we will introduce the model, 
discretize it on the lattice (section B), and derive the equations of motion. In section 3 
we introduce an adapted version of the "Male" and "Female" fermion fields [2] required 
to generate the fermion correlators with c-number fields. In section 4 we describe the 
numerical setup and the results, and we conclude in section 5. 
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2. The Axial-U(l)-Higgs-fermion model in 1+1 dimensions 

We will consider the 1+1 dimensional Abelian-Higgs model, coupled axially to fermions. 
The action reads in the continuum: 

S = Sh + Sa + Sf, (2.1) 

in terms of the components 

Sh^-J d^x [D^(t>'^D^^ct> + A(0t0 - v''/2f] , (2.2) 

Sa^- j d'x -^F^^F^^ (2.3) 

Sf^- j d^x {d^ + zA^75) ^ + Gi, {(P*Pl + <PPr) ^] , (2.4) 

and with the definitions 

D^(P ^ d^(P - iA^(P, F^, ^ d^A, - d,A^, Pfl,i = ^ (1 ± 75) . (2.5) 
The action is invariant under gauge transformations of the form 

^ exp{iq^{x)^5)^, (j) ^ exp{-i^{x))(j), ^ - ^^^(a;), (2.6) 
if we take q — 1/2. 

In lattice simulations it is more convenient to work with vector gauge symmetry, 
rather than axial, and so noting that the left and right chiral components have opposite 
charge, it is therefore natural to charge-conjugate one of them [1], 

^«=(^^C)^, ^^ = -(CV;,)^, ^L=^L, ^L=^L, (2.7) 

where C is the charge-conjugation matrix given in Appendix A. Upon doing this, the 
action in the new variables (but omitting the primes) reads 



J d^x [V^t'' {d^ - iqA^) i/j + \G'il^^C^(j>*'il; - \G^l^C(t>i)'^] . (2.8) 



It is no longer axially coupled, but the Yukawa interaction has become a Majorana 
term, and the gauge symmetry has become vector-like 

ip ^e-KY>{-iqi{x))^, ^ exp(-i^(x))(/), A^ ^ A^,- d^,i{x), (2.9) 
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with the continuum equations of motion being 

D^D''(f) - 2A((/)*0 - v'^/2)(f) - ^t/j^Ct/j = 0, (2.10) 

7'^L>^V' + = 0, (2.11) 

d.F^'' + e\jl^^ + = 0. (2.12) 

We have introduced the gauge currents 

if^) = iqh'^, if;) = i(<t>D^c{>* - <f>*D^<f>). (2.13) 

There is one further symmetry of this system, the one that this work is principally 
interested in, and it is the global U(l) symmetry, — >■ exp{—iuj^^)ip. This symmetry 
has an associated current 

= i^j^^Y^^, (2.14) 

which is precisely the fermion current in the original^ theory, and classically conserved 
if one naively applies the equations of motion. Quantum mechanically, however, it is 
the subject of an anomaly 

« = ^^'"'F,^'^ = 9,C^^ = ^e'^'C., (2.15) 

and this allows us to relate the total fermion number, Q{t) = J dx j^, to the Chern- 
Simons number, C{t) = J dx = f dxAi{x), through 

Q{tf)-Q{ti) = C{tf)-C{U). (2.16) 

There is one further number that is worth mentioning, the winding number of 
the Higgs field. When the Higgs field is away from zero, it takes values on a circle 

parametrized by its phase 9, (f){x) = \(f){x)\e''^^'^\ Using this phase we may define a 
Higgs winding number, describing the number of times the field winds around this 
circle on a given spatial section, 

Nw = ^ I dxdi0{x). (2.17) 
In J 

In a vacuum state we know that the covariant derivative of the Higgs field vanishes, 
and that its modulus is constant, in which case we have that d^d = A^, leading to the 
sum of the Higgs winding and Chern-Simons numbers vanishing in the vacuum. 

For the numerical work, we discretize the Abelian-Higgs-fermion model on a 1+1 
dimensional lattice of size L — aiUx at the level of the action, and derive lattice 
equations of motion as described in Appendix B. 

^Non-charge conjugated. 
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3. Bosons and fermions 

We are interested in the time-evolution of this system, and we will adopt the approach 
of [1], where the dynamics of bosonic and fermionic degrees of freedom are treated dif- 
ferently. The gauge and scalar fields are evolved using the classical equations of motion 
described in the previous section. Classical dynamics is an excellent approximation to 
the quantum dynamics for processes dominated by infrared physics and for fields with 
large occupation numbers. The fermions are treated completely quantum-mechanically, 
in the sense of solving the quantum equation of motion (2.11) in the classical bosonic 
background, in terms of field operators. Since the fermions are bi-linear in the action, 
the equation of motion is linear, and the field can in all generality be expanded in terms 
of a set of mode functions and time-independent creation- annihilation operators (see 
below) . 

This leaves the question of the back-reaction of the fermions on the classical bosonic 
fields. Following [1] again, we interpret the fermionic terms in the gauge and scalar 
equations of motion as expectation values of the corresponding operators, evaluated 
in some state encoded in the expectation values of the creation-annihilation operators. 
These states are time-independent, and amount to specifying an initial condition. The 
time-evolution is in the mode functions only. 

We then take one step further by representing these creation-annihilation operators 
by a set of random numbers, thereby generating an ensemble of fcrmion field-realisation 
[2]. These can each be evolved in the same bosonic background, and the field expec- 
tation values are then replaced by simple averages over the ensemble. The point is to 
note that the number of field realisations (Ng) in the ensemble can be much smaller 
than the number of mode functions (n^), and the statistical approach can therefore be 
much cheaper in terms of computational effort. 

3.1 Boson initialisation 

We will consider two setups for the bosonic fields. The first (in section 4.1) is to by hand 
set the gauge-Higgs evolution to be a sequence of sphaleron transitions, thus forcing 
the Chern-Simons number to change (as in [1]). The fermion fields evolve dynamically 
in the background of these handmade sphalerons. We will use this setup to test the 
ability of the ensemble to capture the anomaly, and to find out how large the ensemble 
needs to be. 

When considering the non-perturbative field dynamics (in sections 4.2 and 4.3), we 
instead initiahse the bosonic fields by setting Aij,{x, t — 0) — 0, 9o0(x, t — 0) — and 
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introducing random noise for the scalar field, 0(x, t = 0) 




(3.1) 



1 2 

in terms of random numbers 0^.' , with the correlator 



2 



(3.2) 



The gauge field momenta OqAi (x, t = 0) are found by solving the Gauss constraint 
(2.12) with fermion sources. 

As described in [21, 22] this initialisation represents^ an initial quantum vacuum 
before Higgs symmetry breaking, Vini = Xv'^cffcf). In the subsequent evolution, momen- 
tum modes /c^ < A^;^ will grow exponentially, and from some time on they can be 
described using classical dynamics. The fermions do not grow, and are still treated 
quantum mechanically. 

The amount of growth of the scalar modes is determined by the (in 1+1 dimensions 
dimensionless) parameter v. This can be seen in various ways. The growth lasts until 
backreaction from self-interactions kick in, i.e. when 0^ ~ v'^. For a given mode, we 
have 



where initially, rik = 0. Classical dynamics is a good approximation once the mode has 
grown so much that + 1/2 ^ 1/2. Hence large v allows for classicality. 

Another way of phrasing this is to note that once ~ the scalar-gauge interaction 
and the effect of the scalar on the fermions goes as ev and Gv, respectively, whereas 
back-reaction of fermions on bosons is e and G. Hence for large fermion effects are 
relatively smaller (the fields have relatively smaller amphtude). 

In the following, we will employ = 64 and v — d). Since only modes with k"^ < Xv"^ 
are unstable, only they will be classical, and these are therefore the only bosonic modes 
we initialise. 

^In fact, we should also initialise the momenta docp with random numbers for the identification with 
the quantum vacuum to be completely correct. Setting dQcj) to zero initially makes the initial total 
charge on the lattice vanish, a requirement for consistency of Gauss law. To achieve this is cumbersome, 
but possible, when initialising both field and momenta. For our purposes here, initialising only the 
field variables will suffice. 
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3.2 Fermion mode expansion 

Now we need to know how to set up the initial conditions for the fermion field, and 
for this we will be using the usual mode expansion. There is a slight complication, 
however, due to the fact that the fermion equation of motion is not linear in -0, but 
involves both -0 and -0* (2.11). This leads to the real and imaginary components having 
different equations of motion, particularly when the lattice Wilson term is included, 
and so it is convenient to write the Dirac spinor as a combination of two Majorana 
spinors, which in our conventions (Appendix A) just means breaking into real and 
imaginary parts. 

V^=i=[^i-i*2], (3.4) 

and it is these components that we write as a mode expansion 

*(^, x)= f^T^ IbkUke''-^ + biVkC-H , (3.5) 
J Zn ZUk L J 

in terms of the constant spinors U and V (given in Appendix A) and a set of creation- 
annihilation operators b^. We then note that the fields are canonically normal- 
ized, and that their conjugate momenta are ^"^12, so the canonical anti-commutations 
relations are 

{*a(t,^),*^(i,^')} = H^-^)^aP, (3.6) 

which may be achieved by imposing 

{^bk,bl,^ = {2n)2ukS{k-k'). (3.7) 

In the equations of motion for the bosonic fields we require the quantum expectation 
value of fermion bilinears, and so we follow [2] in constructing the two-point functions 

D>^{x,y) = (|*«(a;)*^(y)|), D<^{x,y) = -(|*^(y)*«(a;)|), (3.8) 
Da^ix,y) = i [D>^{x,y) + D<^{x,y)] , (3.9) 

leading to 

Dapix, y) = 1 V ^-L [UkaVk^e""-^^-'^ - VkaUkpe-'^-^"""^] . (3.10) 

Z ^ — ' ZTT ZUOk 

where we take 6fc|) = 0. We note that although the fields are real, the two-point function 
is imaginary, D^^{x,y) = —Dap^x^y). The observation of [2] is that we can construct 
a bi-linear of classical spinor fields, for which the ensemble average two-point function 
matches (3.10). This allows us to simulate the quantum backreaction of fermion fields 
using ensemble averages of classical spinor fields; this is what we shall now do. 
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3.3 Fermion ensemble, Male and Female 

If we were to simply evolve an ensemble of fermions, where we draw the initial conditions 
of each realization from a sample with the appropriate distribution and then take the 
ensemble average (\t'(a;)\t'(y)) to mimic the quantum two-point function, we cannot 
reproduce (3.10). However, if one introduces two "genders" of fermions, male and 
female, and writes their mode expansion as 

1 r Hk 1 

*m(i) = ^ y [%U,e"- + 4V,e-' '] , (3.11) 

then we find that taking 

(rjkv;) = i27r)2ujk5{k-p), {rjkVp) = 0, (3.13) 

leads to 

^{^Ma'^'F^) = IJ2 [UkaVkpe''-^'-"^ - VM^e-''-^'-y^] (3.14) 

k 

^ D^p{x,y), (3.15) 

so we now have an explicit way of replacing quantum averages, (l-'^l), with ensemble 
averages, {X) . This leads to us evolving 

Di^D'y - 2A(0> - wV2)0 - !^(^^'^C^^) = 0, (3.16) 

'j'^D^^^^^ + G^V''^'^'* = 0, (3.17) 
d^i&^A^ - dl^A^) + e\f^ + 2l) = 0, (3.18) 

rather than the equations of motion appearing in Appendix B. The fermion gauge- 
current is also modified in this prescription, with the requirement of its conservation 
leading to 

i/-. = 7 [iv;^(x)7;.c/^(x)V'^(x+/i)+iv;^(x+/.)7,c/r(x)v^(x) 

-#^(x + /x)7^[/^*(a;)V"^(x) - ii,\xY, ^Ul{x)i,'\x + /x)] . (3.19) 
Furthermore, we need a representative of the anomalous current, 

Jm,5 = \ [#^(x)7^75C/^(^)V'''(^ + /i)+^VS''(^ + /^)7M75C/;(^)V'^(^) 

-%i^^{x + /x)7^75C/;(x)V^(x) - %i,^{x)^^^^U^{x)i^'^{x + . (3.20) 

Because of cancellation between lattice doublers, this quantity is conserved for vanishing 
Wilson term (ri = 0, see Appendix B), but for r\ — \ the current is anomalous, as we 
will see below. 
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4. Results 



4.1 Hand-made Sphaleron transitions 
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Figure 1: The bosonic fields are evolved 
through a series of sphaleron transitions, 
thereby continuously changing the Chern- 
Simons number. The fermion fields are 
evolved in this background using the equa- 
tions of motion, and fermion number is 
seen to obey the anomaly equation. The 
Higgs winding number follows in steps 
(shown is —N^)- 



Figure 2: For large enough Chern- 
Simons number, the anomaly equation is 
no longer satisfied on a small lattice. In- 
creasing the volume allows a larger range 
of agreement. The three curves coincide 
upon rescaling by the lattice size Ux along 
both axes (insert). 



We first want to check that our approach of replacing mode functions by a random 
ensemble still leads to correct dynamics and that the anomaly equation is satisfied, as 
in [1]. We also want to determine how large the ensemble needs to be to get statistically 
reliable results for the anomaly and the dynamics. 

An elegant way of doing this is to by hand set the gauge-Higgs field evolution to 
be a series of sphaleron transitions, thereby continuously changing the Chern-Simons 
number in a controlled way. The explicit expression for the bosonic fields can be found 
in [1]. The important point is that sphaleron transitions take place at half-integer 
values of t/to, and the fields are in vacuum at integer values. We choose the timescale 
to so that the transitions are slow enough that the fermions, which are evolved using the 
equations of motion in the sphaleron-vacuum background, do not lag too much behind, 
TTiAto = 4. From the point of view of the fermion, the evolution is almost adiabatic, 
and no additional spurious particle creation takes place. Only the particles associated 
with the anomaly contribute. At the sphaleron configuration, the Higgs field length 
vanishes, and Higgs winding changes discontinuously from one integer to the next. In 
the vacuum, A^^ = —Ncs, and we will always plot —N^. 
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Fig. 1 shows the evolution of Chern-Simons number A^^cs, Higgs winding number 
— A\v and the fermion number N{. The parameters used were m^L = evL = 25.6, 
G = 0, Ng = 90, Tlx = 128, V = 2, timestep dt = 0.05. The agreement between 
Chern-Simons number and fermion number is remarkably precise, even for such a small 
ensemble. 

As was pointed out in [1] , fermion number is periodic on a finite lattice with period 
2nx, and so for large Ncs the agreement will fail. Fig. 2 shows fermion number for very 
large Chern-Simons number at different values of the lattice size (volume is fixed 
evL = 6.4). These show the lattice behaviour, and can indeed be rescaled by (along 
both axes) to end up on top of each other (inset). This is exactly as in [1], and means 
that sufficiently large lattices can accommodate any Chern-Simons number. Notice 
that the ensemble is still Ng = 90. 
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Figure 3: Convergence of the fermion 
number as the ensemble is enlarged. Here 
10 to 2430 realisations, for a small lattice 
Ux = 32. 
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Figure 4: Convergence of the fermion 
number as the ensemble is enlarged. Here 
10 to 2430 realisations, for a large lattice 
Ux = 128. 



Although the small ensemble very convincingly reproduces the anomaly when 
looked at by eye, it is only prudent to investigate the statistical precision. This is 
shown in Figs. 3 and 4 for ensembles of 10, 30, 90, 270, 810 and 2430 random re- 
alisations, respectively. The left-hand plot is on a ra^ = 32 lattice, and we see that 
the agreement is always fairly good, even for Nq = 10. Looking closer (inset), we do 
see that the curves converge, and in fact converge to a value slightly off Ncs- This is 
the finite volume effect as described before. In the right-hand plot with Ux = 128, 
this discrepancy is gone and increasing the ensemble, fermion number converges to the 
Chern-Simons number value. We conclude that convergence in Ng is achieved at the 
few-percent level for A^^, = C(IOOO). 

As reported in [1], including the Yukawa coupling G makes the lattice artefacts 
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Figure 5: Fermion number for a small 
lattice rix = 32, with increasing values 
of the Yukawa coupling G. Chern-Simons 
number is shown for comparison. Lattice 
artifacts are larger for non-zero G. 
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Figure 6: Fermion number for fixed G = 
0.1 and volume evL, but increasing Hx, 
decreasing lattice spacing. Also for finite 
Yukawa coupling, lattice artifacts can be 
removed by increasing the number of lat- 
tice points or increasing the volume. 



stronger. Fig. 5 shows the fermion number in the hand-made sphaleron background 
for Ux = 32, evL = 6.4 with varying G. The anomaly holds until A'cs — 5, after 
which the finite size effects kick in, stronger with increasing G. However, increasing 
the lattice size again ameliorates the situation, as shown in Fig. 6, where G/e = 0.1 is 
kept constant, and the lattice discretization is made finer (constant volume evL = 6.4 
and increasing^). 

From a practical point of view, we would 
like to be able to run baryogenesis simula- 
tions for timescales rriAt = evt = C(IOO), 
on a large enough lattice to fit in the ap- 
propriate physics m^L = evL ^ 1, while 
having the anomaly correctly reproduced for 
realistic values of the Chern-Simons number, 
say Ncs = C(10). And we also need to in- 
clude a non-zero Yukawa couphng, at least 
for physics around the electroweak transi- 
tion. The question is whether we can find 
a combination of evL, n^, G and A'"^ that 
can accommodate this. 

Fig. 7 shows a run on a much larger lat- 




Figure 7: Sphaleron transitions and 
fermion number for larger volume evL = 
51.2, at G/e = 0.1. Larger statistics cures 
the discrepancy. 



tice, evL = 51.2, 



2, with Yukawa cou- 



^ Increasing Ux with increasing physical volume has the same effect. 



- 12 - 



pling G/e = 0.1,ina range of A^cs = — 10. We first note that the larger volume makes 
the anomaly agree less well than for the runs in Fig. 6, which did reasonably well until 
Ncs = 5. This is the case both for n^. = 256 and for = 512, with half the lattice 
spacing (not shown). However, this is just adding up of statistical fluctuations, and 
can be compensated for by increasing the ensemble size. At A'^j = 810 the agreement is 
again convincingly reproduced. The finite volume effect is not apparent at these lattice 
sizes. 

4.2 Non-equilibrium dynamics 




Figure 8: The Higgs field, Chern-Simons 
number, winding number and fermion 
number during a tachyonic transition. 
Chern-Simons number and fermion num- 
ber are indistinguishable. 
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Figure 9: Convergence in Nq of Chern- 
Simons and fermion number for v = 64. 
Convergence is excellent until m^t ~ 100, 
after which the effects of statistical fluctu- 
ations in the fermions have accumulated 
enough to make a difference. 



In a fast-quench symmetry breaking transition, Higgs field modes with k'^ < 
will be unstable ("tachyonic") and grow exponentially. This drives the gauge field 
to also grow until non-linear backreaction begins to dominate, stop the growth and 
eventually leads to thermalisation. In the presence of CP-violation, such a transition 
may lead to a baryon asymmetry (see for instance [13]). 

For our choice of initial conditions, the initial gauge field is driven by the fermion 
ensemble fluctuations and the initial scalar fleld. Our goal is that for a given scalar fleld 
conflguration, the evolution should be independent of Nq, so that the statistics reliably 
reproduce the fermion state. We need a large enough ensemble to have statistical 
fluctuations under control. We will set 

A 1 

6'^ 4 



256, mA 



ev 



0.2, ruH = V2Aw2, 



(4.1) 
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and vary Nq. 

Fig. 8 shows the Higgs field (black line), Higgs winding number (green), Chern- 
Simons number (red) and fermion number (blue) in such a transition. G — 0, v — 6A 
and Nq — 2430. The Higgs field "falls off the hill" as expected, and performs oscil- 
lation around its finite temperature minimum. Meanwhile, the Chern-Simons number 
grows and oscillates near the integer-value Higgs winding number. The anomaly is so 
well obeyed that fermion number is essentially indistinguishable from Chern-Simons 
number. The Higgs winding bounces around in the beginning, an effect of the Higgs 
field length being small, and the winding number therefore ill-defined. But once sym- 
metry breaking gets going, winding number is stable, integer and consistent with the 
Chern-Simons number. 

We illustrate the convergence of the dynamics with increasing Nq in Fig. 9, where 
tachyonic transitions are performed for v — 64 for different sizes of the ensemble. As 
expected, we see convergence in Nq, but also that the required ensemble is C(IOOO), to 
get agreement at this value of v and for these times. In fact, the Chern-Simons number 
is very sensitive to fiuctuations in the fermion source. This is at least partly because in 
1 + 1 dimensions, a U(l) gauge field only has one dynamical degree of freedom (i.e. up 
to gauge transformations), which is precisely the Chern-Simons number. This means 
that in the sea of fermion degrees of freedom, the single degree-of-freedom gauge field 
can easily be bounced around. These issues are specific for 1-1-1 dimensions, and we 
will proceed with v — 6A and Nq — 2430, for which convergence is under control at 
least for niAt < 100, and quahtatively correct for m^i < 150. This will suffice for the 
present work, but can be improved depending on the level of precision required. 

Since the Yukawa coupling is absent, the fermions are massless throughout. Also, 
the fermion and boson total charges are individually conserved at the level of O{10~^^), 
and Gauss law is conserved at a (relative) level of O{10^^). 

For non-zero Yukawa coupling the fermions acquire masses as the Higgs transition 
proceeds (in addition to the gauge and Higgs fields). As we saw in section 4.1, the 
Yukawa coupling introduces stronger lattice artefact, but these could be cured by using 
larger ensembles. In Fig. 10 and 11 we show the evolution and convergence in Nq of 
a simulation with G/e — 0.1, v = 64. We see that although there is a clear effect of 
non-zero G, convergence still holds by increasing the ensemble to a few thousand, and 
fermion number (dashed lines) follows Chern-Simons number fairly well. 

4.3 Application: Cold Electroweaik; Bciryogenesis in 1-|-1 dimensions 

In the minimal model of electroweak baryogenesis, CP-violation is provided through 
the CKM fermion mass matrix. For hot baryogenesis, this effect is much too small 
to account for the asymmetry and a separate source of CP-violation is required. The 



- 14 - 




Figure 10: Chern-Simons number, Higgs Figure 11: Convergence in Nq of Chern- 
winding number, fermion number and Simons and fermion number for f = 64 at 
Higgs field in a tactiyonic transition at G/e = 0.1. At finite Yukawa coupling, we 
G/e = 0.1, Nq = 2430, v = 64. need a somewhat larger ensemble to reach 

convergence, here Nq = 2430 — 7290. 
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Figure 12: The ensemble averages 
of Higgs field, Chern-Simons number, 
fermion number and Higgs winding num- 
ber over 64+64* scalar realisations, with 
V = M and k = 0.04. 
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Figure 13: A zoom-in of Fig. 12. The 

asymmetry is driven by the oscillation of 
the Higgs field, through the C(P) vio- 
lating force term. Winding number can 
only change when the (local) Higgs field 
is small. 



situation is less clear for "cold" baryogenesis (see for instance [13]). For illustration, we 
will postpone this issue, and simply introduce C(P) violation^ in our 1+1 dimensional 



"'From the point of view of the present model, we actually break C and P separately, while CP is 
conserved. This is the analogue of requiring CP violation in 3+1 dimensions. 
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model through a bosonic term in the action (exactly as in [21]), 

S^S- J d'^x^ (t>*(t> e^,F^\ (4.2) 
which amounts to an addition to the bosonic equations of motion (2.12), (2.10) of 

dlA,^... + ^^do\<j>\\ (4.3) 
9o^0 = . . . + ^^0^10, (4.4) 
d,d^A^ = ...-'^d,\cp\\ (4.5) 

The conservation of Gauss law and the anomaly and the convergence in Nq is unaltered 
by this addition, and the C(P) violation is not obvious from a given random realisation 
of the bosonic fields. We now need to also average over an ensemble of bosonic real- 
isations, each with a separate ensemble of fermion fields. This represents a quantum 
initial state of the Higgs fields coupled to fermions initially in the vacuum. 

Fig. 12 shows the scalar-ensemble averaged obscrvablcs, Higgs field, Chern-Simons 
number, fermion number and winding number, aX k = 0.04, w = 64 and Nq = 2430. 
We average over a set of 64 random realisations plus the corresponding C(P)-conjugate 
configurations. This makes the ensemble explicitly C(P) symmetric, and the asymmetry 
will be identically zero ior k — 0. This procedure is similar to the one employed in 
[24, 25]. We see that an asymmetry is indeed generated in Chern-Simons, winding and 
fermion numbers as the transition proceeds. The anomaly is very well obeyed until 
times m^t ~ 100, where fermion number goes a little low. We checked that this is 
indeed due to the configurations with relatively large A^cs = 8 — 10, for which the lattice 
artefact makes a small deviation. 

From Fig. 13, a blow-up at early times, we see that because Higgs winding can 
only take place in the presence of a local Higgs field zero, the asymmetry is created 
when the average Higgs field is low in its oscillation. When it is high, winding number 
is essentially constant. Also, since the C(P)-violating force is proportional to the 
gauge field picks up speed in-between Higgs extrema. The net effect is a "pumping" 
behaviour, of the baryon asymmetry as Higgs symmetry breaking proceeds. As the 
Higgs approaches a uniform vev, and oscillations damp out, the asymmetry creation 
gradually stops. This is very similar to the case in the 3+1 dimensional SU(2)-Higgs 
model [12, 24, 25], where the gauge field is much more complicated dynamically than 
the model considered here. 

One important question is to what extent fermions can be ignored dynamically, 
compared to the bosonic fields. If not, many simulations of baryogenesis may need 
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Figure 14: The evolution of average 
Chern-Simons number, with (black) and 
without (red) fermion back-reaction. Here 
shown for v = 64 (small effect) and v = 8 
(some effect). 
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Figure 15: Dependence of the asym- 
metry on C(P)-violation k, measured at 
rriAt = 150. The behaviour is linear for 
small K. Note that the asymmetry is zero 
at «; = by construction. 



corrections from the fermion backreaction. In Fig. 14 we show the average Chern- 
Simons and fermion number with and without fermion backreaction, for f = 64 and 
V = 8, respectively. Clearly, for v = 64, fermion backreaction can be mostly ignored, 
and the fermions only serve as spectator fields, encoding the fermion asymmetry. For 
this case, one might as well just do bosonic simulations, and infer fermion number from 
the anomaly equation a posteriori. For v = 8 however, the fermions begin to influence 
the evolution of the gauge field, even when the bosonic fields are subject to a tachyonic 
instability and therefore grow large. For the Standard Model in 3+1 dimensions, there 
is no f -ambiguity, and it will be crucial, to what extent back-reaction is important. In 
particular, CP-violation itself is a backreaction effect, which will dynamically generate 
effective terms similar to the bosonic C(P)-violation used here. 

Finally, to illustrate the type of calculations that are possible, we show how the 
asymmetry depends on k, (Fig. 15). For small enough k, the dependence is nicely 
linear. Also remember that because our scalar ensemble is explicitly C-symmetric, 



5. Conclusion 

By combining the methods of [1] and [2], we have demonstrated how to do first-principle 
numerical simulations of bosonic scalar-gauge systems with quantum fermions, in a nu- 
merically efficient manner. Although there is no gain in numerical effort in the specific 
1+1 dimensional toy model considered here (compared to [1]), in the physically relevant 



-17- 



case of 3+1 dimensions, we expect a significant decrease in the required computing time. 
As an example, including fermions in the simulations of [12], an ensemble of Nq — 2430 
should be compared to the 90^ (= 729, 000) mode functions otherwise required, a gain- 
factor of 300. Either way, fermions are numerically challenging, but for the setup of 
[12], there would then be no need for the resource-consuming CP- violating term. Ex- 
cept for the scalar ensemble averages in section 4.3, all the simulations presented here 
were done on a normal desktop computer in less than 24 hours in total. 

We found that Gauss' law, fermion back-reaction as well as the baryon anomaly 
are well reproduced in terms of a statistical ensemble of (9(1000) fermion field reali- 
sations. As described in [2] implementation of the fermion correlators, including the 
anti-commutativity of the fermionic operators requires a doubhng of the fields into 
"male" and "female" (not to be confused with the standard lattice fermion doublers), 
adapted to the system at hand. In addition, the usual lattice doubler problem has to 
be addressed; in the present case we found that Wilson fermions in space and a small 
timestep was sufficient to keep the doublers sufficiently decoupled that they stayed un- 
excited for the timescales required here. Failure to do this leads to an exact cancelling 
out of the baryon anomaly. 

The method requires careful consideration of the interplay between lattice size, 
ensemble size and the size of couplings. In particular, the Yukawa coupling introduces 
additional lattice artefacts, which have to be compensated for, and we have demon- 
strated how to do this. 

The upshot is that fermions are included completely in the dynamics, since they are 
bi-linear in the action, and so at least in cases where gauge fields are dominated by large 
particle numbers and long wavelength, this approach provides a very reliable description 
of the full field dynamics. The obvious application of the method is (electroweak) 
baryogenesis, where baryon number violating processes are classical in nature, whereas 
the CP-violation^ and the actual baryon number are carried by the fermion degrees of 
freedom. This applies both to "hot" and "cold" baryogenesis. 

Simulations of the early stages of the heavy-ion collisions are also within the scope 
of the work presented here, since it involves very large (boosted) gluons fields coupled 
to (sea and valence) quarks. The valence quarks source the gauge field, which then 
evolves and may in turn source the emission of fermions. With the method presented 
here, fermions may be included in the dynamics completely. 

The obvious next step is to implement ensemble fermions in 3-1-1 dimensions, cou- 
pled to SU(2)-Higgs bosonic fields as in the Standard Model, where the gauge-fermion 
interaction is chiral rather than axial. Including the Standard Model CP-violation via 

^At least in the Standard Model. 



-18- 



the CKM matrix will require all three fermion generations, and represents a signifi- 
cant numerical challenge; with the method described here, this numerical effort can be 
reduced by one or even two orders of magnitude. This is work in progress. 



A. Conventions 

Wc use the metric signature (—,+), and for the Dirac algebra, we employ the Weyl- 
Majorana representation 

{7^7'^} = 2<^ VS^^V'V, 75 = -7°7i, (A.l) 

with explicitly 

-'-'^^^ — v^(j:),(a.) 

C=(°7), r^ = -CrC-\ (A.3) 

B. Lattice equations 

On the lattice, we define the gauge link variables (a^, /i = 0, 1 are the lattice spacings), 

U^{x) = exp{-ia^A^{x)), U^{x) = exp(-iga^^^(x)). (B.l) 

We are using the non-compact formulation of the gauge action, with the basic gauge 
field variable. We define the derivatives 

d.A^ix) = - {A,{x + /x) - A,{x)) , D^ = hD, + D'^] , (B.2) 

D^(f> = —[U^{x)cl>{x + //) - = -Ulix- pl)cI>{x - //)], (B.3) 

D,^{x) = -[U^^{x)^{x + /x) - ^{x)], D'^^ix) = -[^l;{x) - C/^t(^ _ ^)^(^ _ ^)], 

(B.4) 

We will deal with the spatial fermion doublers by including a Wilson term 

WiiP = -|riaiL>iL>V. (B.5) 

The lattice action then becomes, 

Si^^^^Sh + Sa + Sf, (B.6) 
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with 

Sh = Y1 [^00^^00 - ^10^^10 - A(0t0 - ^2/2)2] , (B.7) 



x,t 



= E ^ (^0^1 (^) - 9iAo{x)f , (B.8) 

X 

Sf^-J2 «oai [ip {h'"i^^^ + D'^) + W)iP + \G(t)*VC^^P - \G(t)^pC^P^] . (B.9) 

x,t 

This immediately gives the lattice equations of motion 

Df^D'^cj) - 2X{(f)*(f) - vy2)(j) - ^VCilj = 0, (B.IO) 

^^b,i^ - ^^D,D[^ + G(t>r = 0, (B.ll) 

d,{&^A^ - &^A^) + e^{j) + 3l + 3w) = 0, (B.12) 

and the currents, 

jfe,^ = 2(0D^0*-0*D^</.), (B.13) 

= f [^l^{xh,U^^{x)^^;{x + pt) + i^{x + ii)^^m;{x)^{x)\ , (B.14) 

fw = 0, (B.15) 

J^^ = [^"^ - , (B.16) 

and so we have Gauss' law from 
The chiral current is 

X — — 

3n,5 = 2 ['^i^hf^^5U^,{x)ipix + + ip{x + ii)-f^-f5U*{x)ip{x)] , (B.18) 
= 0. (B.19) 

Chern-Simons number and fcrmion number is trivially adapted to the lattice, and 
following [23], we write = |0(a;)|e*^^^\ and then define the integer lattice winding 
number as 

a; 
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Figure 16: Tachyonic transition run at 
V = 6A, K = 0, = 270 with different 
values of tlie Wilson coefficient ri. With- 
out the Wilson term, the spatial doublers 
cancel out the anomaly, ri = 1 seems a 
good choice. 
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Figure 17: The same tachyonic tran- 
sition as on the left, now with different 
values for the timestep dt. Even at the 
largest timestep, the timelike doublers are 
not sufficiently excited to cancel out the 
anomaly. We use the smalest timestep 
shown here, dt = 0.05. 



C. Fermion doublers 



Doublers contribute to the anomaly with the opposite sign to the non-doubler modes, 
and the anomaly will then average out if we do not remove the doublers from the 
dynamics. Fig. 16 shows the anomaly in a tachyonic transition for "naive" fermions 
ri = and with a Wilson term ri > 0. The anomaly disappears for the naive fermions, 
when doublers are allowed to get excited. 

By adding a Wilson term in space, but not time, the fermion equation will still 
lead to temporal doublers, so where we thought we were evolving a single Fermi-field 
ip, we are actually evolving two Fermi-fields, which we call and To see this we 
define 



if t is even, 
if t is odd. 



(c.i; 



Now, if the lattice time derivative is evaluated on an even t slice, then it actually only 
uses fields evaluated on the preceeding odd t slice, and the following odd t slice, in 
which case one finds the equation of motion (B.ll) becomes 







airi 



DiD[^+{x) + G(f){x)ilj+*{x) 



(C.2) 



-7i 



airi 



YD^^-{x) - -^DiD[tlj-{x) + G<j){x)tp-*{x) 
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Similarly, on odd t slices we have 




(C.3) 



+7i 



YD^ip-{x) - ^lDiD[i^-{x) + G<P{x)ij~*{x) 



showing that both ■0+ and satisfy the fermion equation of motion, and that we are 
actually evolving two fermi degrees of freedom. 

In Fig. 17 we sec a set of runs where we vary the timestep. We only initialise the 
non-doubler modes, and we see that for all timcstcps that give a stable integration of 
the equations of motion, the time-like doublers stay un-excited. In all simulations in 
the main paper, we use dt = 0.05, the smallest time-step presented here, and we see no 
doubler effects. 



D. Spinors 

We now construct the basis spinors required in the mode expansion of the fermion 
operators, which we do by setting Ai — and taking = so that for 

ruf — Gv/\/2. For the positive frequency modes 



*i = C/i,ike^'=", A;° = a;i>0, (D.l) 



we are then led to 



U,.= (-'^'), (D.2) 

1 Ti 

Sk — — sin(ai/c), = — [1 — cos(ai/c)], (D.3) 
Mi^mk + mf, c^i = + V(Mi)2 + (sfe)2, (D.4) 

Similarly, the negative frequency solutions are found by 

*i = ^l,fee-^^•^ = a;i > 0, (D.5) 

but because the field is Majorana, and therefore real, we immediately have V — U*. In 
order to calculate the two-point function we need the following identity, 

Ui,kUi,k - Ml - ir~k^, (D.6) 

where — {ui, s^). To get the mode functions for ^f^, simply make the replacement 
ruf -nif. 
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